/*******************************************************************************
	   Author: 	Neil Himwich
	   Project: ETS
	   Purpose: Create PLANT-MONTH, PLANT-WEEK, PLANT-DAY, PLANT-PERIOD KG
	   EMISSIONS, using the MOC IMPUTATION RULE (ie the penalization system
	   used in market). Save these as rule "M" for "MOC"
	   
	   Note: This adapted the code found in this directory from other
	   imputation rules, and from the archive for the MOC imputation rule

	   Date created: 4 March 2024
	   Version: STATA 18 MP
	   
	   Edited by: Noah Sobel-Lewin 
	   Date: 22 August 2024
   ****************************************************************************/


set more off
clear all
pause on

*-------------------------------------------------------------------------------

use "$EMISSIONS_DATA_IN/CEMS_2018-12-31 to 2021-04-03_StackDay_Balanced_HistCFs_Master337s_3sample2trunc.dta", clear
drop if date < date("2019jan16", "YMD")
format date %tddd-Mon-YY
sort date composite_id
tab date
** 242,303 Observations. 337 Stacks. 719 Days from 16-Apr-19 to 03-Apr-21.

** Generate a week variable
drop week
gen int week = floor((date-td(31dec2018))/7)+1
la var week "Week"

** Generate a month variable which goes from the 16th of each month to the 15th of the next month. 
** This definition of month aligns with the implementation of ETS market in its initial Stages. 
** Nick Cox's code:
drop month16
gen month16 = date if day(date) == 16
// 	replace month16 = mdy(1, 16, 2019) in 1
replace month16 = month16[_n-1] if missing(month16)
la var month16 "Month, defined as 16th of this month to 15th of next month"
gen month16_sif = month16
format month16 %tddd-Mon-YY
la var month16_sif "[SIF] Month-year. For convenience."
order month16_sif, a(month16)
sort composite_id date

** Generate Dummy for Interregnum period of market
** Covid lockdown -- 23 March 2020 to 11 October 2021. 
** Mock-III started on 12 October 2021.
gen D_interregnum_week = week >= 65 & week <= 93
replace D_interregnum_week =1 if week == 99 | week == 100
la var D_interregnum_week ///
	"[Week-level] =1 if Interregnum due to Covid-19 lockdown, zero otherwise."

drop extinct_with_assignment removed_from_sample comment

** Remove plants that have no observation from 16-Apr-19 to 03-Apr-21.
bysort gpcb_id: egen pct_report = mean(!missing(pm_mass_rule0)) if ///
	month16 >= date("2019jul16", "YMD") & D_interregnum_week == 0
sort composite_id date
bysort gpcb_id: egen pct_report_1 = mean(pct_report)
drop if pct_report_1 == 0
drop pct_report pct_report_1
sort composite_id date  
tab date
** 222,890 Observations. 310 Stacks. 719 Days from 16-Apr-19 to 03-Apr-21.

gen period = 0
** Mock-I
replace period = 1 if date >= td(15jul2019) 
** Mock-II
replace period = 2 if date >= td(13aug2019)
** Comp-I
replace period = 3 if date >= td(16sep2019)
** Comp-II
replace period = 4 if date >= td(16oct2019)
** Comp-III
replace period = 5 if date >= td(16nov2019)
** Comp-IV
replace period = 6 if date >= td(01jan2020)
** Comp-V
replace period = 7 if date >= td(01feb2020)
** Comp-VI
replace period = 8 if date >= td(01mar2020)
** Interregnum-I
replace period = 9 if date >= td(22mar2020)
** Mock-III
replace period = 10 if date >= td(12oct2020)
** Interregnum-II
replace period = 11 if date >= td(12nov2020)
** Comp-VII
replace period = 12 if date >= td(01dec2020)
** Comp-VIII
replace period = 13 if date >= td(01jan2021)
** Comp-IX
replace period = 14 if date >= td(01feb2021)
** Comp-X
replace period = 15 if date >= td(01mar2021)

drop if date >=  td(01apr2021)

*tab period if composite_id == "GJSRT100740111"
gen temp = 1
bysort period: egen period_length = sum(temp) 
replace period_length = period_length / 310 
drop temp

********************************************************************************
********************* [1] CALCULATE ACTUAL STACK DAY MASS **********************
********************************************************************************

forvalues i=0(1)9{

	** Replace tot_report_hrs to zero 
	** if the stack was not sending a valid caliberated stack-day kg/hr that day
	gen tot_report_hrs_rule`i' = tot_report_hrs
	replace tot_report_hrs_rule`i' = 0 if D_vctransmitting_rule`i' == 0
	label var tot_report_hrs_rule`i' ///
		"[RE-LABELED] Stack-day VALID CALIBRATED Reporting Hours from CEMS - rule `i'"

	** Calculate DAILY actual mass emissions using CEMA Day Reporting Hours
	gen stack_daily_mass_actual_rule`i' = pm_mass_rule`i' * tot_report_hrs_rule`i'
	la var stack_daily_mass_actual_rule`i' ///
		"Stack Day Actual Mass Emissions (kg) from CEMS Reported Hours - rule `i'"

}

********************************************************************************
***************** [2] Weekly NON-REPORTING HOURS FOR IMPUTATION *****************
********************************************************************************

forvalues i=0(1)9{

	replace tot_report_hrs_rule`i' = 0 if tot_report_hrs_rule`i' == .
	bysort composite_id week: egen tot_weekly_report_hrs_rule`i' = sum(tot_report_hrs_rule`i')
	bysort composite_id week: gen weekly_DA_rule`i' = (tot_weekly_report_hrs_rule`i'/168)*100
	gen non_report_hrs_rule`i' = 168-tot_weekly_report_hrs_rule`i'
	by composite_id week: egen mean_weekly_report_hrs_rule`i' = mean(tot_report_hrs_rule`i')
	** Calculate DAILY non-reporting hours for imputation
	gen daily_nonreport_hrs_rule`i' = 24 - tot_report_hrs_rule`i'
	order daily_nonreport_hrs_rule`i', a(tot_report_hrs_rule`i')
	la var daily_nonreport_hrs_rule`i' ///
		"Stack Day NON-Reporting Hours - CEMS - rule `i'"

}

********************************************************************************
********* [3] For MOC Imp: **********
********************************************************************************
*encode composite_id, gen(composite_id_enc)
gen composite_id_enc = composite_id
forvalues i=0(1)9{

// Note that the percentiles and means are calculated at the week level, which 
// is the appropriate level to match the imputation method used in the market.

	** Calculate stack's own week mean kg/hr
	bysort composite_id week: egen imp_mean_mass_rule`i' = mean(pm_mass_rule`i')
	label var imp_mean_mass_rule`i' "Stack Week Mean kg/hr - rule `i'"
	
	** Calculate stack's own week p25 kg/hr
	bysort composite_id week: egen imp_p25_mass_rule`i' = pctile(pm_mass_rule`i'), p(25)
	label var imp_p25_mass_rule`i' "Stack Week P25 kg/hr - rule `i'"
	
	** Calculate stack's own week p50 kg/hr
	bysort composite_id week: egen imp_p50_mass_rule`i' = pctile(pm_mass_rule`i'), p(50)
	label var imp_p50_mass_rule`i' "Stack Week P50 kg/hr - rule `i'"
	
	** Calculate stack's own week p75 kg/hr
	bysort composite_id week: egen imp_p75_mass_rule`i' = pctile(pm_mass_rule`i'), p(75)
	label var imp_p75_mass_rule`i' "Stack Week P75 kg/hr - rule `i'"
	
	** Calculate stack's own week p90 kg/hr
	bysort composite_id week: egen imp_p90_mass_rule`i' = pctile(pm_mass_rule`i'), p(90)
	label var imp_p90_mass_rule`i' "Stack Week P90 kg/hr - rule `i'"
	
	tempfile helper_p90
	save `helper_p90', replace
	
	** Calculate stack's own 90 day p90 kg/hr leading up to this compliance period
	tsset composite_id_enc date
	mvsumm pm_mass_rule`i', stat(p90) win(90) gen(threemnthp90_mass_rule`i') force
	bysort composite_id period: keep if _n == 1
	keep composite_id period threemnthp90_mass_rule`i'
	merge 1:m composite_id period using `helper_p90', nogen
	
	** If all last three months missing then assign it the 8.08 we use when no data reported 
	replace threemnthp90_mass_rule`i' = 8.08 if missing(threemnthp90_mass_rule`i')
	
}


********************************************************************************
******************** [4] CALCULATE WEEKLY/DAILY STACK IMPUTED MASS ********************
********************************************************************************

local tot_weekly_hrs = 168
local week_5pc = 0.05*168 // for 100-95% DA bin
local week_15pc = 0.15*168 // for 95-80% DA bin
local week_30pc = 0.30*168 // for 80-50% DA bin
local week_1pc = .01*168 // for less than 1% DA

local week_95pc = .95*168
local week_80pc = .8*168
local week_50pc = .5*168

forvalues i=0(1)9{	

	***Categorize industry into imputation bin based on weekly data availability/non-reporting hours
	***Create a bin variable called cat = 0,1,2,3,4 with 5 if 100% weekly DA, and 0 if <1% weekly DA.
	egen cat_rule`i' = cut(tot_weekly_report_hrs_rule`i'), at(0, `week_1pc', `week_50pc', `week_80pc', `week_95pc', `tot_weekly_hrs', 169) icodes
	
	***Calculate imputed mass quantity based on graded imputation bin, using cat
	// If report full week
	gen stack_weekly_mass_imp_rule_mkt_`i' = 0 if cat_rule`i'== 5
	// If have DA between 100 and 95
	replace stack_weekly_mass_imp_rule_mkt_`i' = non_report_hrs_rule`i'*imp_mean_mass_rule`i' if cat_rule`i' == 4
	// If have DA between 95 and 80
	replace stack_weekly_mass_imp_rule_mkt_`i' = `week_5pc'*imp_mean_mass_rule`i' + ((non_report_hrs_rule`i'-`week_5pc')*imp_p75_mass_rule`i') if cat_rule`i' == 3
	// If have DA between 80 and 50
	replace stack_weekly_mass_imp_rule_mkt_`i' = `week_5pc'*imp_mean_mass_rule`i' + (`week_15pc'*imp_p75_mass_rule`i') + ((non_report_hrs_rule`i'-`week_5pc'-`week_15pc')*imp_p90_mass_rule`i') if cat_rule`i'== 2
	// If have DA between 50 and 1
	replace stack_weekly_mass_imp_rule_mkt_`i' = `week_5pc'*imp_mean_mass_rule`i' + (`week_15pc'*imp_p75_mass_rule`i') + (`week_30pc'*imp_p90_mass_rule`i') + ((non_report_hrs_rule`i'-`week_5pc'-`week_15pc'-`week_30pc')*threemnthp90_mass_rule`i') if cat_rule`i'== 1
	// If have DA less than 1
	replace stack_weekly_mass_imp_rule_mkt_`i' = 168*8.08 if cat_rule`i'== 0 | missing(tot_weekly_report_hrs_rule`i')
	
	drop imp_p*_mass_rule`i'
	foreach pct in 50 55 60 65 70 75 80 85 90 95 {
		bysort composite_id: egen imp_p`pct'_mass_rule`i' = pctile(pm_mass_rule`i') if ///
		(date >= td(16Jul2019)) & (D_interregnum_week == 0), p(`pct')
		
		gen stack_weekly_mass_imp_rule_p`pct'_`i' = non_report_hrs_rule`i'*imp_p`pct'_mass_rule`i'
		
	}
		
	
*** B) Calculate PER DAY imputed emissions
	foreach imprule in mkt p50 p55 p60 p65 p70 p75 p80 p85 p90 p95 {
			***TRANSLATE WEEKLY IMPUTED QUANTITY TO DAILY IMP QUANTITY BY DIVIDING IT BY 7. THIS IS MY WORKAROUND FOR DAILY APPLICATION OF WEEKLY MOC RULE.
			***This is problematic, but StackDayPMMassRuleM.dta is nowhere used in the analysis.
			gen stack_daily_mass_imp_rule_`imprule'_`i' = stack_weekly_mass_imp_rule_`imprule'_`i'/7
			la var stack_daily_mass_imp_rule_`imprule'_`i' ///
				"Stack-day Imputed Mass Emissions (kg) from CEMS Non-Reported Hours - rule`i'"

			** Correct for week of Diwali 2019 when imputation should be zero.
			replace stack_daily_mass_imp_rule_`imprule'_`i' = 0 if ///
				date >= date("2019oct26", "YMD") & ///
				date <= date("2019nov03", "YMD") & ///
				stack_daily_mass_imp_rule_`imprule'_`i' != .

	}

}

********************************************************************************
******************* [5] CALCULATE DAILY STACK VALIDATED MASS *******************
********************************************************************************

forvalues i=0(1)9{	

	foreach imprule in mkt p50 p55 p60 p65 p70 p75 p80 p85 p90 p95 {
		
		** Calculate validated DAY emissions per stack
		gen stack_daily_mass_val_rule_`imprule'_`i' = stack_daily_mass_actual_rule`i' + stack_daily_mass_imp_rule_`imprule'_`i'
		replace stack_daily_mass_val_rule_`imprule'_`i' = stack_daily_mass_imp_rule_`imprule'_`i' ///
			if stack_daily_mass_val_rule_`imprule'_`i' == .
		la var stack_daily_mass_val_rule_`imprule'_`i' ///
			"Stack Day Validated Mass Emissions (kg) (=Stack-day Actual + Stack-day Imp) - rule `i'"
	}

}

/*******************************************************************************
	   NOW READY TO MAKE THE DEPENDENT VARIABLE FOR REGRESSIONS IN PART 6 
	   (6. TREATMENT EFFECT ANALYSIS). THIS VAR IS PLANT MONTH MASS (KG).

	   STEPS TO CREATE THIS VARIABLE WILL BE:

	   6.1) CALCULATE STACK WEEK MASS (KG), AND THEN PLANT WEEK MASS (KG).
	   FIRST GENERATE THESE VARIABLES IN THE STACK-DAY LEVEL DATA SET.
	   AND SAVE IT FOR FUTURE USE IF NEEDED ("_StackDayLevel.dta").

	   6.2) THEN, COLLAPSE THE STACK-DAY LEVEL DATA SET TO KEEP ONE DATA POINT 
	   PER PLANT PER WEEK ("_PlantWeekLevel.dta"). 

	   CONTAINS THE DEPENDENT VARIABLE: PLANT-WEEK MASS (KG) EMISSIONS 
	   USING IMPUTATION RULE 0a.
   ****************************************************************************/

preserve
keep date composite_id gpcb_id treatmentstatus stack_daily_mass_actual_rule0 stack_daily_mass_imp_rule*0 stack_daily_mass_val_rule*0
order date composite_id gpcb_id treatmentstatus stack_daily_mass_actual_rule0 stack_daily_mass_imp_rule*0 stack_daily_mass_val_rule*0
rename stack_daily_mass_actual_rule*0 stack_daily_mass_actual_*M
rename stack_daily_mass_imp_rule*0 stack_daily_mass_imp_*M
rename stack_daily_mass_val_rule*0 stack_daily_mass_val_*M 
sort date gpcb_id composite_id
save "$IMPUTATION_DATA_OUT/StackDayPMMassRuleM.dta", replace

restore


/*******************************************************************************
	   NOW READY TO MAKE THE DEPENDENT VARIABLE FOR REGRESSIONS IN PART 6 
	   (6. TREATMENT EFFECT ANALYSIS). THIS VAR IS PLANT MONTH MASS (KG).

	   STEPS TO CREATE THIS VARIABLE WILL BE:

	   6.1) CALCULATE STACK MONTH MASS (KG), AND THEN PLANT MONTH MASS (KG).
	   FIRST GENERATE THESE VARIABLES IN THE STACK-DAY LEVEL DATA SET 
	   AND SAVE IT FOR FUTURE USE IF NEEDED ("_StackDayLevel.dta").

	   6.2) THEN, COLLAPSE THE STACK-DAY LEVEL DATA SET TO KEEP ONE DATA POINT 
	   PER PLANT PER MONTH ("_PlantMonthLevel.dta"). 

	   CONTAINS THE DEPENDENT VARIABLE: PLANT-MONTH MASS (KG) EMISSIONS 
	   USING IMPUTATION RULE 1a.
   ****************************************************************************/

preserve
   
********************************************************************************
********************* [6] GENERATE STACK-MONTH MASS (KG) *********************
********************************************************************************

forvalues i=0(1)9{	
	
	foreach imprule in mkt p50 p55 p60 p65 p70 p75 p80 p85 p90 p95 {

		** Calculates MONTHLY validated emissions PER STACK
		bysort composite_id month16: egen stack_month_mass_val_rule_`imprule'_`i' = sum(stack_daily_mass_val_rule_`imprule'_`i')
		bysort composite_id month16: egen miss_stack_month_mass_val_`imprule' = max(missing(stack_daily_mass_val_rule_`imprule'_`i'))
		replace stack_month_mass_val_rule_`imprule'_`i' = . if miss_stack_month_mass_val_`imprule' == 1
		la var stack_month_mass_val_rule_`imprule'_`i' ///
			"Stack-month Validated Mass Emissions (kg)"
		drop miss_stack_month_mass_val_`imprule'
	
	}

}	

** Keep one monthly data point for each STACK
bysort composite_id month16: gen dump=_n
bysort composite_id month16: keep if _n==_N 	
keep composite_id industry_id gpcb_id month16 treatmentstatus stack_month_mass_val* 
** 7,440 observations. 310 stacks. 24 months from April 2019 to March 2021.

******************** [6b] COLLAPSE TO PLANT-LEVEL DATA SET ********************

/**	AS PER NICK'S INSTRUCTION ON WEEKLY CALL 01 DEC 2020: FOR PLANTS WITH >1 STACK, 
	WE WILL ONLY CALCULATE THEIR PLANT MASS (KG) SUM IF ALL STACKS ARE CONNECTED. 
	IF STACKS ARE PARTIALLY CONNECTED (i.e., ONLY 1 OF 3 STACKS CONNECTED), 
	THE ENTIRE PLANT SUM WILL BE MISSING. 	**/

forvalues i=0(1)9{
	
	foreach imprule in mkt p50 p55 p60 p65 p70 p75 p80 p85 p90 p95 {

	** Calculate MONTHLY validated emissions PER PLANT
	bysort gpcb_id month16: egen ind_month_mass_val_rule_`imprule'_`i' = sum(stack_month_mass_val_rule_`imprule'_`i')
	bysort gpcb_id month16: egen miss_ind_month_mass_val_`imprule' = max(missing(stack_month_mass_val_rule_`imprule'_`i'))
	replace ind_month_mass_val_rule_`imprule'_`i' = . if miss_ind_month_mass_val_`imprule' == 1
	la var stack_month_mass_val_rule_`imprule'_`i' ///
		"Plant-month Validated Mass Emissions (kg) - rule `i'"
	drop miss_ind_month_mass_val_`imprule'
	
	}
	
}

** Keep one monthly data point for each PLANT
bysort gpcb_id month16: gen dump=_n
bysort gpcb_id month16: keep if _n==_N
keep industry_id gpcb_id month16 treatmentstatus ind_month_mass_val* 
** 7,008 observations. 292 plants. 24 months from April 2019 to March 2021.

save "$IMPUTATION_DATA_OUT/PlantMonthPMMassRuleM.dta", replace

restore

preserve

********************************************************************************
********************* [7] GENERATE STACK-WEEK MASS (KG) **********************
********************************************************************************

forvalues i=0(1)9{

	foreach imprule in mkt p50 p55 p60 p65 p70 p75 p80 p85 p90 p95 {
		
		** Calculate WEEKLY validated emissions PER STACK
		bysort composite_id week: egen stack_week_mass_val_rule_`imprule'_`i' = sum(stack_daily_mass_val_rule_`imprule'_`i')
		bysort composite_id week: egen miss_stack_week_mass_val_`imprule' = max(missing(stack_daily_mass_val_rule_`imprule'_`i'))
		replace stack_week_mass_val_rule_`imprule'_`i' = . if miss_stack_week_mass_val_`imprule' == 1
		la var stack_week_mass_val_rule_`imprule'_`i' ///
			"Stack-week Validated Mass Emissions (kg) - rule `i'"
		drop miss_stack_week_mass_val_`imprule'
		
	}

}

** Keep one monthly data point for each STACK
bysort composite_id week: gen dump=_n
bysort composite_id week: keep if _n==_N 	
keep composite_id industry_id gpcb_id week treatmentstatus stack_week_mass_val*
** 31,930 Observations. 310 Stacks. 103 weeks from April 2019 to March 2021.

******************** [7b] COLLAPSE TO PLANT-LEVEL DATA SET ********************

forvalues i=0(1)9{

	foreach imprule in mkt p50 p55 p60 p65 p70 p75 p80 p85 p90 p95 {
		
		** Calculate WEEKLY validated emissions PER PLANT
		bysort gpcb_id week: egen ind_week_mass_val_rule_`imprule'_`i' = sum(stack_week_mass_val_rule_`imprule'_`i')
		bysort gpcb_id week: egen miss_ind_week_mass_val_`imprule' = max(missing(stack_week_mass_val_rule_`imprule'_`i'))
		replace ind_week_mass_val_rule_`imprule'_`i' = . if miss_ind_week_mass_val_`imprule' == 1
		la var stack_week_mass_val_rule_`imprule'_`i' ///
			"Plant-week Validated Mass Emissions (kg) - rule `i'"
		drop miss_ind_week_mass_val_`imprule'
		
	}

}

** Keep one monthly data point for each PLANT
bysort gpcb_id week: gen dump=_n
bysort gpcb_id week: keep if _n==_N
keep industry_id gpcb_id week treatmentstatus ind_week_mass_val* 
** 30,076 Observations. 292 Plants. 103 weeks from April 2019 to March 2021.

save "$IMPUTATION_DATA_OUT/PlantWeekPMMassRuleM.dta", replace

restore

preserve 

********************************************************************************
********************* [8] GENERATE STACK-PERIOD MASS (KG) ********************
********************************************************************************

forvalues i=0(1)9{
	
	foreach imprule in mkt p50 p55 p60 p65 p70 p75 p80 p85 p90 p95 {

		** Calculate MONTHLY validated emissions PER STACK
		bysort composite_id period: egen stack_period_mass_val_rule_`imprule'_`i' = sum(stack_daily_mass_val_rule_`imprule'_`i')
		bysort composite_id period: egen miss_stack_period_mass_val_`imprule' = max(missing(stack_daily_mass_val_rule_`imprule'_`i'))
		replace stack_period_mass_val_rule_`imprule'_`i' = . if miss_stack_period_mass_val_`imprule' == 1
		la var stack_period_mass_val_rule_`imprule'_`i' ///
			"Stack-period Validated Mass Emissions (kg) - rule `i'"
		drop miss_stack_period_mass_val_`imprule'
	
	}

}

** Keep one monthly data point for each STACK
bysort composite_id period: gen dump=_n
bysort composite_id period: keep if _n==_N 	
keep composite_id industry_id gpcb_id period period_length treatmentstatus stack_period_mass_val* 
** 3,410 observations. 310 stacks. 11 periods.

******************** [8b] COLLAPSE TO PLANT-LEVEL DATA SET ********************

/**	AS PER NICK'S INSTRUCTION ON WEEKLY CALL 01 DEC 2020: FOR PLANTS WITH >1 STACK, 
	   WE WILL ONLY CALCULATE THEIR PLANT MASS (KG) SUM IF ALL STACKS ARE CONNECTED. 
	   IF STACKS ARE PARTIALLY CONNECTED (i.e., ONLY 1 OF 3 STACKS CONNECTED), 
   THE ENTIRE PLANT SUM WILL BE MISSING. 	**/

forvalues i=0(1)9{

	foreach imprule in mkt p50 p55 p60 p65 p70 p75 p80 p85 p90 p95 {
		
		** Calculate MONTHLY validated emissions PER PLANT
		bysort gpcb_id period: egen ind_period_mass_val_rule_`imprule'_`i' = sum(stack_period_mass_val_rule_`imprule'_`i')
		bysort gpcb_id period: egen miss_ind_period_mass_val_`imprule' = max(missing(stack_period_mass_val_rule_`imprule'_`i'))
		replace ind_period_mass_val_rule_`imprule'_`i' = . if miss_ind_period_mass_val_`imprule' == 1
		la var ind_period_mass_val_rule_`imprule'_`i' ///
			"Plant-month Validated Mass Emissions (kg) - rule `i'"
		drop miss_ind_period_mass_val_`imprule'
		
	}

}

gen temp = 1
bysort gpcb_id: egen num_stack = sum(temp)
replace num_stack = num_stack / 16
drop temp

** Keep one monthly data point for each PLANT
bysort gpcb_id period: gen dump=_n
bysort gpcb_id period: keep if _n==_N
keep industry_id gpcb_id period period_length treatmentstatus num_stack ind_period_mass_val* 
keep if period > 0
** 2,920 observations. 292 plants. 16 periods.

save "$IMPUTATION_DATA_OUT/PlantPeriodPMMassRuleM.dta", replace
